Random Matrix Theory approach to Mesoscopic Fluctuations of Heat Current 
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We consider an ensemble of fully connected networks of N oscillators coupled harmonically with 
random springs and show, using Random Matrix Theory considerations, that both the average 
phonon heat current and its variance are scale-invariant and take universal values in the large 
iV-limit. These anomalous mesoscopic fluctuations is the hallmark of strong correlations between 
normal modes. 



PACS numbers: 44.10.+i, 66.10.cd, 64.60.ae, 63.20.-e 

Introduction- The study of heat conduction by 
phonons in disordered or chaotic structures have at- 
tracted recently considerable interest [TH3]- A central 
issue of these investigations is the dependence of the av- 
erage heat current J on the system size N. A naive 
expectation is that disorder or phonon-phonon interac- 
tions scatters normal modes and induces a diffusive en- 
ergy transport that leads to a normal heat conduction 
described by Fourier's law which states that J ~ TV -1 . 
Many studies [MS], however, find that in low dimensional 
chains J scales as J ~ N~ a , where a is usually differ- 
ent from one. In fact, experiments on heat conduction 
in nanotubes and graphene flakes have reported observa- 
tions of such anomalous behavior [TJHdT] . 

However, many real stuctures such as biological sys- 
tems [12] and artificial networks in thin-film transistors 
and nanosensors |13j are not simple one-dimensional or 
two-dimensional lattices. Rather they are characterized 
by a complex connectivity that can be easily designed 
and realized in the laboratory [TH - Hrj] . Therefore, not 
only is it a fundamental demand for the development of 
statistical physics to understand normal and anomalous 
heat conduction in complex networks of coupled oscilla- 
tors, but it is also of great interest from the technological 
point of view, since the achievements of modern nano- 
fabrication technology allow us to access and utilize such 
structures with sizes in the range of a few nanometers up 
to few hundred nanometers. 

The complexity of coherent wave interferences in such 
networks calls for a statistical treatment of any of their 
transport characteristics. This way of thinking has been 
adopted already in classical wave and quantum transport 
theories associated with mesoscopic chaotic or disordered 
systems and resulted in a plethora of exciting results 
[TTRT^] like the weak and strong Anderson localization, 
the universal conductance fluctuations (UCF) etc. The 
statistical approach led also to the revival of Random 
Matrix Theory (RMT) |20J, as a major theoretical tool 
for the analysis of transport characteristics of complex 
systems. RMT has found applications in many areas 



of physics ranging from nuclear, atomic and molecular 
physics to mathematical physics (for a review see [21]). 
Consequently a variety of RMT ensembles have been in- 
troduced [22] , extending the original work of Wigner be- 
yond the traditional Gaussian ensembles, helping to un- 
derstand phenomena like the Quantum Hall effect |23j . 
Anderson localization and Metal-to-Insulator transition 
[21] . The success of RMT was such that in recent days it 
has become almost a dogma that this theory captures the 
universal properties of complex disordered and chaotic 
systems [i~7Hl9] [21] . It is thus surprising, that the study 
of fluctuations and the use of RMT as a concrete tool for 
their analysis were not brough up in any of the previous 
studies of heat transport. 

In this Letter we address heat transport and the asso- 
ciated sample-to-sample mesoscopic fluctuations of com- 
plex networks of N equal masses connected with one 
another via random harmonic springs. The force ma- 
trix that describes the dynamics of the system is real 
symmetric and consists of random elements (spring con- 
stants). We find that the statistical description of heat 
transport can be effectively described by an ensemble of 
Random Matrices with diagonal elements that fluctuate 
with a variance N times larger than the corresponding 
variance of the off-diagonal elements. Using RMT consid- 
erations we show that both the average heat current (J) 
and its variance (A J) 2 are scale-invariant and get a uni- 
versal value in the large- N limit. These anomalous meso- 
scopic fluctuations is the hallmark of strong correlations 
between normal modes of the system. For moderate size 
networks, with random springs taken from a distribution 
with variance a 2 < 1/N, we find that the heat transport 
is sensitive to the boundary conditions imposed on the 
two end-sites which are coupled to the thermal baths. 
In particular, for fixed boundary conditions the current 
is completely dominated by a pair of surface modes, for 
which only the end sites oscillate with appreciable ampli- 
tude. We hope that our analysis will motivate the use of 
RMT models and provide new insight in the mesoscopic 
fluctuations of heat transport. 



Fully Connected Harmonic Networks — We consider a 
network of N harmonic oscillators of equal masses m = 
mQ. The system is described by the Hamiltonian [25 



n = ] i p t m- 1 p+] i q t ^q 



(1) 



where Q T = (qi,q 2 ,--- ,gjv), pT = (Pi,P2,--- ,Pn) and 
q n ,Pn are respectively the individual oscillator displace- 
ments and momenta. The mass matrix is M nm = 5 nm mo, 
and $ is the force matrix that contains also information 
about the boundary conditions (b.c). For a fully con- 
nected network of coupled oscillators with free b.c. $ 
takes the form $ nm = ( J2i k n i)5 nm - k nm where k nm are 
the spring coupling constants. These spring constants 
knm are chosen to be symmetric {k nm = k mn ) and uni- 
formly distributed according to k nm £ [— ^ + 1, ^ + l] 
where the disorder strength parameter W has to be 
smaller than 2 in order to ensure that k nm > 0. In the 
case of fixed b.c. $ has to be modified by considering the 
coupling of the first and last oscillator to hard walls i.e. 

^nm = ®nm + (koiSnlSml + fciVTV+l^niV^miVb 

Next, we want to study the non-equilibrium steady 
state (NESS) of this network driven by a pair of Langevin 
reservoirs set at temperatures T® and T® (we assume 
T B > T B ), and coupled to the first n = 1 and last n = N 
masses with a constant coupling strength 7. The corre- 
sponding equations of motion that describe also the cou- 
pling to the bath are q n = dH/dp n , p n = —dH/dq n + 

(-JPn/m + ^2jT®( n ) (<5„i + 5 nN ), where £„(*) is 
delta-correlated white noise ( n (t)( n '(t') = S nn 'S(t — t'). 



The NESS current is evaluated as J 



-(Tf-Ji) 



— (TV — T B ) where the temperature of the n-th oscil- 
lator is defined as T n = p n /mQ. The notation 7TT which 
will be implicetly assumed from now on, indicates the 
thermal statistical average. 

For weak coupling 7 it was shown in Ref. p] [2] that 



j = J2 jM; J [M) = Cv 
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where /„ = \ipn \ 2 and ipn indicates the n— th com- 
ponent of the /x-th normal mode of the Hamiltonian Eq. 
IT) and the coefficient C = ^(Tf - T 2 B ) [H]. In Eq. 
2), the /i-th addendum J^' is naturally interpreted as 
the contribution of the /x-th mode to the total heat flux 
J. As intuitively expected, J^ is larger for modes that 
have larger amplitudes at the boundaries and couple thus 
more strongly with the reservoirs. Thus the analysis of 
heat flux J reduces to the study of the normal modes of 
Hamiltonian % given by Eq. (fl). 

Random Matrix Theory Formulation - We separate 
out the random component of the spring constants and 



re-write them as k r . 



1 + W nm where now W n 
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FIG. 1: (Color online) Free boundary conditions, (a) A typi- 
cal distribution of rescaled heat flux J = J/ Co for a network 
of N — 10 3 oscillators and disorder strength W = 1. (b) The 
rescaled average heat current (J) (open symbols) and vari- 
ance (AJ) 2 = (AJ) 2 /Cq (filled symbols) versus a. Various 
system sizes iV (indicated in the figure) have been used. The 
dashed lines in (b) are the results of D-RMT ensemble with 
strongly fluctuating diagonal elements. Both in (a) and (b) 
we have used Eq. (pf in order to evaluate the heat current J. 



a constant matrix A and a random part R as 

<£ = A + R where A = ATI - U; R = D-W (3) 

where 1 is the N x N unit matrix, U is a matrix whose 
all elements are equal to unity i.e. U nm = 1, I? is a 
diagonal matrix with D nn — —^2i^ n W n i, and W is a 
random matrix defined below. The above decomposi- 
tion allow us to distinguish the various contributions. 
The matrix W can be treated as a "standard" RMT 
ensemble (note though that it has zero diagonal ele- 
ments) . It is convenient to rewrite it as W = uWq where 
Wq is a RM with elements having unit variance where 
a 2 = (AW nrn ) 2 = W 2 /12. The diagonal matrix D has 
Gaussian distributed random elements with (D nn ) = 
and variance (AD nn ) 2 = (N — l)a 2 . 

The constant matrix A can be diagonalized exactly. It 
has (a) one eigenvalue Wo = with a corresponding eigen- 
vector (l/\/F)(l, 1, 1 • ■ ■ , 1) T and (b) N - 1 degenerate 
eigenvalues w M = N (/i = 1,2, ■•■TV— 1). Now consider 
adding to A the random matrix W, i.e. we neglect for 
the moment D and consider 



$' = A + aW 



(4) 



w w 
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] . The force matrix <3& can be decomposed into 



Already for an arbitrary small a, the N — 1 degeneracy 
will be removed and the corresponding eigenvectors will 
be those of a (N — 1) x (N — 1) random matrix. The 
(N — l)-time degenerate level is broadened into a band 
of width ~ o-y/N. The perturbation theory applies for 
ay/N < N i.e. for er < \/N. However even for larger a 



the RMT still applies because then we can simply neglect 
the matrix A in Eq. W. In short, for small a we have 
an RMT for (N — l)-rank matrices (the contribution to 
current of the level with ujq « can be neglected, in 
comparison to the (N — l)-levels), whereas for large a we 
have an RMT for iV-rank matrices. Thus, in the large 
N limit we treat Eq. Q as an ensemble of N x N GOE 
matrices [17] . 

Normalization requires that (if ) = (1$ ) = (I n ) = 

1/N. Defining a rescaled variable Xn = In /(I n ), we 
can rewrite Eq. (pi) as 

Co 
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(5) 



According to the standard RMT, and omitting the 
mode label /^, the joint probability distribution of 
the rescaled eigenmode intensities X n is a prod- 
uct of two Porter-Thomas distributions P(X\,Xn) = 
(l/2%)(l/y/X 1 X N ) exp[-(Xj + X N )/2]. Assuming fur- 
ther that the various z^' -terms appearing in Eq. (15J) are 
statistically independent we get 



(J) = -Co; (AJ) 2 



8N 



r 2 



(6) 



Comparison of these theoretical predictions with a direct 
numerical evaluation of the mean and the variance of 
heat current J via Eq. (pi) (see Fig. IIId) leads us to 
conclude that standard RMT considerations describe well 
the scaling of the average current but not the variance. 
To obtain the correct description of the variance, it is 
necessary to treat the full force matrix, as given in Eq. 
([3]). Below we show that the matrix D induces strong 
correlations between different z^'s, thus, invalidating 
the assumption which led to Eq. (pjl) for the variance 

D-RMT ensemble with strongly fluctuating diagonal el- 
ements - We now consider the ensemble of matrices 
given by Eq. (I3J). Again for large N, the matrix A has 
no effect, so it is enough to understand the eigenvectors 
of the random matrix R. The eigenvalues of D are of 
order \D nn \ ~ avN, so that they occupy a band of or- 
der ay/N and are separated by a typical energy interval 
aVN/N = cr/yN. The same is true for the eigenvalues 
of the matrix W. In this sense D and W are "of the same 
strength" and neither can be treated as perturbation to 
the other. However, the qualitative understanding of the 
eigenvectors of the combined matrix R is along the fol- 
lowing lines: The eigenvectors of D are localized on the 
individual sites i.e. the fi— th eigenvector is ipn — S nu . 
The matrix W mixes these eigenvectors, so that eigen- 
vectors of R are spread over all sites and resembles those 
of a standard RMT. Therefore (J) is qualitatively not 
different from the standard RMT result of Eq. (f6j) . The 
only difference is that the coefficient 1/4 now assumes 
the numerical value ss 0.19 (see Fig. [TV)) . 



As far as the variance (A J) 2 is concerned we get results 
that are qualitatively different from the standard RMT 
result of Eq. (pi). It turns out that in this case each 
eigenvector of R "remembers" the set (Z?n,--- ,Dnn) 
of the eigenvalues of D so that correlations between dif- 
ferent eigenvectors of R are significantly stronger than 
those for the standard RMT. Namely the mode-mode 
correlations between the different ?M's of the matrix R 
are described by [27J 

(« M « M ) - (z (ri )(z W )(l + e) = (z) 2 (l + e) (7) 

where e is a constant. Using Eq. ([7| we calculate the 
variance (AZ) 2 of the random variable Z (see Eq. (Km)): 



(AZ) 2 = N 2 e(z 2 )+0(N). (8) 

Expressing (A J) 2 in terms of Z via Eq. (pM) we get 

(AJ) 2 = C 2 e(z) 2 , (*)«0.19 (9) 

Direct numerical evaluation of the variance (A J) 2 based 
on Eq. (pi) confirms the above theoretical estimates. 
In Fig. [1J3 we show some of our numerical results for 

rescaled variance (AJ) 2 = (AJ) 2 /Cq. The data clearly 
indicate that (AJ) 2 is scale invariant for any disorder 
strength a. Further 1/iV numerical analysis allow us to 
extract the asymptotic value e w 0.075. 

We have also checked that correlations between the 
matrices D and W do not play a role in our arguments. 
Detail numerical analysis indicates that if instead of the 
actual D (i.e. D nn = — J^Wn,) we consider a diago- 
nal random matrix completely independent of W so that 
(AR nm ) 2 = er 2 [l + (N — l)S nm ], we still obtain the same 
behavior for (J) and (A J) 2 (dashed lines in Fig. [Tp) . We 
remark that this kind of ensembles, with strongly fluctu- 
ating diagonal elements, (D-RMT ensembles) have pre- 
viously appeared in the context of mescoscopic physics 
[28]. 

Fixed b.c. - Finally we investigate the effect of b.c. on 
the statistics of heat flux. We consider the other limit- 
ing case of fixed b.c. We assume that the first and the 
last oscillator are coupled to the left and right walls with 
spring constants fcoi = I+W01 and fcjvAf+i = 1 + Wnn+i 
respectively which are taken from the same ensemble 
of random springs as the ones in the bulk of the net- 
work. The random components are then included in 
the matrix elements D\\ and Dnn, respectively. The 
constant matrix A also changes to A Rx — A + C where 



(- a 



Snifimi + 5 n NO~mN- This results in a slight shift 



of the zero mode ujq = of the matrix A together with 
a "deformation" of the (1, 1 • • • , 1) T eigenvector. Contri- 
bution of this level to the total current is of order 1/N, 
and it is disregarded below. 

In addition two new levels emerge from the N — 1 de- 
generate subspace of the matrix A. The first one has the 
highest energy u>n-i — ./V+l with a corresponding eigen- 
mode V (Ar_1) = (l/\/2)(l, 0, • ■ • , -1) T . This is an exact 
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FIG. 2: (Color online) Fixed boundary conditions, (a) The 
rescaled average heat current (J) = (J) /Co versus disorder 
strength for various system sizes N. In the main panel we 
scale the a;-axis as aN 3 ' 2 while in the inset (RMT domain) we 
scale it as aN 1 ' 2 . (b) The same as in subfigure (a) but now 
for the rescaled variance (A J) 2 = (A J) 2 /Cq. The various 
symbols correspond to different system sizes N as indicated 
in Fig. 111). The dashed lines are the predictions of D-RMT 
(see also FigHTs) while the solid lines represent the theoretical 
predictions of Eq. |11|12|. 



eigenvalue and eigenvector of A . The second level is 
slightly lower than N + 1 (approximately by 2/N) and 
its eigenvector is symmetric i.e. ujn-2 ~ N+1 — 2/N with 
iP N ~ 2 » (l/\/2)(l - 1/N, -2/N, ■ • • , -2/N, 1 - l/N) T . 
Below we refer to these states as "surface" modes. 

It turns out that for a network described by the con- 
stant force matrix A , most of the current is carried 
by the two surface modes. Using Eq. (pi) we find that 

j(N-l) = j(N-2) = 1 Cq _ At the same tim6j thc rc _ 

maining N — 3 degenerate modes does not contribute to 
the current (in the large N— limit) . Since any of these 
eigenvectors ipv*) has to be orthogonal to both ■0( JV_1 ) 



(a0 



and ip( N 2 ) we get that i[)\ 

These two constrains are satisfied simultaneously only if 



^ and ^ 



-# 



n 



U') 



^ ] = 



for any /i = 1, • • • , N — 3. Thus the total 



heat current is 



J 



N 



:Cn 



(10) 



The above result will still hold as long as the random 
matrix R does not destroy the pair of states tj}^ N ~ x ' and 
ij)( N ~ 2 \ As a increases we observe a coupling of the two 
states towards a linear combination i.e. (l/\/2) \tp^ N ~ r ' ± 
rp( N ~ 2 ']. The origin of this re-organization is traced to 
the matrix D which in the {ip^ N ~ 1 \il)^ N ~' 2 ^} subspace, 
would produce a pair of eigenvalues separated by a dis- 
tance of order o~\fN. This has to be compared to the 
separation of order 1/N between the surface mode eigen- 
values ujn-i and wjv-2 of the matrix A Rx . When a 



reaches a value a c ~ TV" 3 / 2 the two "surface" eigen- 
states are destroyed giving rise to a set of new modes 
that have components (0, • • • ,0, 1) T and (1,0, ••• , 0) T 
i.e. they are localized asymmetrically at the reservoir 
sites. Consequently, the average current will drop to 
approximately a zero value. As the disorder continues 
to increase, the matrix W lifts the degeneracy of the 
N — 3 levels centered around w = N and creates a spec- 
tral band of size Sw ~ cr^/N. For some critical value of 
a = ctrmt ~ 1/yN the bandwidth 6w becomes as broad 
as the gap that separates the degenerate states from the 
surface states. The latter now merge with the continuum 
of states in the band, and the RMT results are recovered. 
For disorder strength such that the dominant contri- 
bution comes only from the two surface states, a quanti- 
tative description of the heat transport can be achieved 
by considering a simple two level system. The two sur- 
face states of the perfect system are described (in the 
site representation) by the 2x2 matrix A^ — — jjU^ 

where U^ has unit elements Unm — 1. This matrix has 
eigenvalues uj\ = 0, L02 — — 2/N and corresponding eigen- 
vectors V (1) = (l/\/2)(l, -1) T and V (2) = (l/\/2)(l, 1) T . 
The two energy levels are separated by an interval 2/N 
where we have set the energy of the highest level (asso- 
ciated to thc the ujn-i level of the original problem) to 
zero. We now add the diagonal matrix D^ with elements 
D^l = VNW nm S nm where W nm G [-W/2,W/2]. The 
total "Hamiltonian" takes the form l>( 2 ) = A^ + D {2 \ 
We can diagonalize exactly this two-dimensional matrix 
and get the corresponding eigenvectors. Using Eq. (pi) we 
obtain J2 = , , ^mrSf^ w ■>? ■ From this we can further 



4+N 3 (W 11 -W 22 ) 

calculate the average and the variance of heat current 
We get 



<J 2 ) = 2C - 



>arctan[f] - log [l + ( 



f) 2 ] 



w = N 3/2 W 
(11) 



while for the variance we get 



(AJ 2 



; arctan [f ] - 8 (log [l + (f ) 2 ] - warctan [f ]f 



2w 4 



r 2 



(12) 

These theoretical predictions are compared in Fig. [2]with 
the numerically evaluated average heat current and vari- 
ance via Eq. (pi) for various system sizes N and disorder 
strength W. Obviously Eqs. ( 11|12 | do not apply for 
ctrmt ^ N~ x / 2 when RMT dominates the transport. 

Conclusions - In conclusion, we have employed RMT 
modeling as a valuable tool for the analysis of mesoscopic 
fluctuations of heat current J in complex (chaotic) net- 
works. For the most basic chaotic system consisting of a 
fully connected network of random springs we have found 
that both the average heat current and its variance are 
scale-invariant. For large AT— limit, these quantities as- 
sume a universal value which is independent of the spe- 
cific boundary conditions. Our analysis indicated that 



the statistical properties of J are affected by the existence 
of correlations between normal modes. For moderate size 
networks with random springs taken from a distribution 
with variance a 1 < 1/iV, the mean and the variance of 
heat current are affected by the existence of two sur- 
face modes emerging in the presence of fixed boundary 
conditions. It would be interesting to investigate the sta- 
tistical properties of heat current for other geometries 
beyond the zero-dimensions, or in the presence of anhar- 
monicities [TJI2S] and establish analogies with mesoscopic 
phenomena observed in the realm of electron transport. 
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Supplementary Material: NESS for the fully 
Connected Network 

In order to establish that the fully connected network 
of harmonic oscillators Eq. (IT]) reaches the NESS, we 
have also performed independent Molecural Dynamics 
(MD) simulations for both free and fixed boundary con- 
ditions. Since these simulations are time consuming we 
confine ourselves to moderate iV-sizes. In Fig. [3] we re- 
post such representative simulations for a case of a fully 
connected network of N = 5 coupled oscillators with 
random springs k nm taken from a uniform distribution 
knm € [1 — W/2; 1 + W/2] and compare these results 
with the ones coming from a direct diagonalization of 
the associated force matrix $ with the use of Eq. (|2| . 

In Fig. [3] open symbols correspond to the average heat 
current and full symbols to its variance evaluated from 
the MD simulations, while the solid lines are the results 
of the diagonalization method that makes use of Eq. (pi) . 
For the MD simulations we have used typically 100 dis- 
order realizations (this has to be compared to the diag- 
onalization method where typically we had more than 
5000 realizations). An additional time average (over the 
last 20 time units) was performed in order to average 
out the oscillations of the chain elements. In order to 
check the convergence of the MD simulations, we have 
compared the flux J for two different times (the time t is 
measured in units of mean inverse frequency) . A conver- 



gence towards the theoretical results of Eq. |2| is evident 
indicating that our system reached a NESS. 
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FIG. 3: (Color online) Molecular Dynamics (MD) simula- 
tions (open and filled symbols) for the case of a network of 
N — 5 fully connected oscillators. The results from the MD 
are compared with the results coming from Eq. d2p. A nice 
agreement, both for the mean heat current (upper) and the 
variance (lower) is observed, indicating that our system can 
reach a NESS. 



